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Abstract. We show that Kahan's discretization of quadratic vector fields is equivalent to 
a Runge-Kutta method. When the vector field is Hamiltonian on either a symplectic vector 
space or a Poisson vector space with constant Poisson structure, the map determined by this 
discretization has a conserved modified Hamiltonian and an invariant measure, a combination 
previously unknown amongst Runge-Kutta methods applied to nonlinear vector fields. This 
produces large classes of integrable rational mappings in two and three dimensions, explaining 
some of the integrable cases that were previously known. 



1. Introduction: Kahan's method for quadratic vector fields 

Consider a system of differential equations arising from a quadratic vector field 

x = f(x) := Q(x) + Bx + c, ieR", (1) 

where Q is an R™-valued quadratic form, B 6 M" xn , and c £ M. n . Consider the numerical 
integration method ii->i' with step size h given by 

= Q(x,x') + ^B(x + x')+c (2) 

where 

Q(x, x') = X - (Q(x + x') - Q(x) - Q{x')) (3) 

is the symmetric bilinear form obtained from the quadratic form Q by polarisation. We call ^ 
the Kahan method. It is symmetric (i.e., self-adjoint), and, crucially, it is only linearly implicit, 
that is, x 1 can be computed by solving a single linear system (because the right hand side of 
([2| is linear in x'). The method ^ was introduced by Kahan in [5] for two examples, a scalar 
Riccati equation and a 2-dimensional Lotka-Volterra system ([8], p. 14) and written down in 
the general form ^ in [5j (see also references therein). 

Because of the different treatment of each term and the unusual treatment of the quadratic 
term, Kahan called ^ an 'unconventional' method. 

The map obtained from applying the Kahan method to various quadratic vector fields / has 
been shown to be completely integrable in a number of cases (see [5], [6], [10] . and references 



Geometric properties of Kahan's method 



2 



therein). In most cases, the conserved quantities depend on the step size h. At present there is 
no single 'integrability mechanism' known which accounts for all integrable cases. 

In this paper we show that the Kahan method is a Runge-Kutta method. As such it shares 
a number of features with all Runge-Kutta methods: it has a B-series, it is affine covariant, and 
it preserves all affine symmetries and all linear integrals of / automatically. As a symmetric 
linear method it preserves all affine reversing symmetries of / automatically, and the B-series of 
its modified vector field contains only even powers of h. 

We then consider the case that / is a Hamiltonian vector field on either a symplectic vector 
space or a Poisson vector space with constant Poisson structure, in any dimension n. We show 
that in this case the Kahan map has a conserved quantity that converges to the Hamiltonian 
of the vector field as h — > 0. It also has a conserved measure which converges to the Euclidean 
measure as h — > 0. These general properties explain some of the integrable cases considered in 

HDJ. 

2. Kahan's method as a Runge-Kutta method 

Proposition 1. The Kahan method coincides with the Runge-Kutta method 

— i/w+v(^)-i/M w 



x — x 
h 

restricted to quadratic vector fields 



Proof. We have 



X r X = Q(x, x') + ^B{x + x') + c 

= \ (Q(x + x') - Q{x) - Q{x')) + \b{x + x') + c 

(5) 

x + x'\ . ,.\ 1 „„x + x' 1 



2 ( 4Q (^) " Q{X) ~ Q{X ' ] ) \ BX + W 



□ 



(Many other Runge-Kutta methods also coincide with the Kahan method when restricted 
to quadratic vector fields. In this paper, we restrict our attention to Q.) As already noted by 
Kahan [9], the Kahan method also coincides with a certain Rosenbrock method on quadratic 
vector fields, for expanding in Taylor series about x gives 



\f(x) + 2(/0r) + \f'{x){x' -x) + \f"{x){x' -x,x'- x)) 
- \ {fix) + f'(x)(x' ~ x) + \f"{x){x' -x,x'- x)) 



(6) 



so 



f{x) + \f{x){x'-x) 

I-p'WY'fix). (7) 
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From the symmetry of the method, or by expanding instead around x', Kahan's method can also 
be written 



(8) 



The B-series of this method is 



°° h k+i 



fe=0 



that is, it contains only tall trees. (For nonquadratic vector fields, the methods Q and Q are 
not necessarily equivalent.) 

The Runge-Kutta method Q has 3 stages and Butcher tableau 
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The modified vector field of the Kahan method applied to quadratic vector fields can be calculated 
using standard methods [4 . Its first few terms are 



/r 



h 4 



f+ 12 ( _2/ "tf> V + f'f'W + 2^0 ( 3/ ' 4/ " 2/ ' 2/ " (/ ' /} 

- 6/7" (/, ff) - 8/"(/> i" 2 /) + 12/" (/, /"(/> /)) + 4/"(/'/, /7)) + • • • • 

A calculation using conjugation by B-series, considering only quadratic vector fields, now yields 
the following result. We omit the details. 

Proposition 2. Kahan's method applied to general quadratic fields has order 2 and is conjugate 
to symplectic up to order 4- It is not conjugate by B-series to a method of order greater than 2 
or conjugate- symplectic by B-series to order higher than 4- 



3. Conservative properties of Kahan's method 

We now consider the conservative properties of the Kahan method in the case of canonical 
Hamiltonian systems x = J _1 V H{x) where H : W L — > R is the Hamiltonian or energy of the 
system. First, note that the method Q is the a = — 1/2 member of the class of Runge-Kutta 
methods 

= af(x) + (1 - 2a)/(^) + af{x'). (9) 

These are all symmetric, A-stable, and second order. Some other members of this family are also 
known to have conservative properties: 

(i) When a = 0, we have the midpoint rule. It is symplectic for canonical Hamiltonian systems. 
Because it is symplectic, it conserves the Euclidean measure. When the Hamiltonian H is 
analytic, the method has a formal invariant H = H + J2T= i h k Hk. When H is quadratic 
(i.e. when / is linear) this series converges to give a conserved quantity of the method. 

(ii) When a = 1/2, we have the trapezoidal rule. It is conjugate to the midpoint rule (the 
conjugacy being half an Euler step), and so it is also conjugate to symplectic and hence 
conserves a measure close to the Euclidean measure, and it also has a formal invariant close 
to H. 
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(iii) When a = 1/6, we have 'Simpson's method' [2], so-called because the right hand side of 
|9| is Simpson's quadrature of A = f(£x + (1 — £)x') d£, appearing in the average vector 
field method ^~E^ = A which conserves the Hamiltonian in canonical Hamiltonian systems. 
Simpson's method preserves quartic Hamiltonians exactly because it coincides with the 
average vector field method in that case. It is not conjugate to symplectic in the sense of 
B-series [I]. 

Proposition 3. Kahan's method has a conserved quantity given by the modified Hamiltonian 
H(x) :=H{x)+ l -hVH{x) T (l- l -hf{x)y 1 f(x) (10) 

for all cubic Hamiltonian systems on symplectic vector spaces and on all Poisson vector spaces 
with constant Poisson structure. The modified Hamiltonian is (i) a rational function of x; (ii) an 
even function ofh; and (iii) given by a convergent series of elementary Hamiltonians containing 
only even-order tall trees. 

Proof. We first consider the homogeneous case, i.e., we let / = KVH(x) where K is an arbitrary 
(not necessarily invertible) constant antisymmetric matrix and H{x) — C(x,x,x) where C is a 
symmetric trilinear form. Note that VH(x) T v — 3C(x,x,v) for all x, v E M. n . For any of the 
methods (J9j|, we have (writing x — (x + x')/2) 

= h (aVH{x) + (1 - 2a)VH{x) + aVH{x')f K T (aVH(x) + (1 - 2a)VH{x) + aVH(x')) 
= h (af(x) + (1 - 2a)/ (5) + af(x') f (aVH(x) + (1 - 2a)VH(x) + aVH{x')) 

= (x' - x) T {aVH(x) + (1 - 2a)VH(x) + aVH(x')) 
3 

= 3aC(x, x, x' — x) + — (1 — 2a)C(x + x', x + x' , x' — x) + 3aC(x' , x', x' — x) 

= | ((2a + l)(C(x', x', x') - C(x, x, x)) + (6a - l)(C(x, x, x') - C(x, x', x'))) . 

The case a = i is Simpson's method, confirming that H (x) is conserved in that case. 



For Kahan's method, a = — ^, and we have from Eqs. (7 



that 



x — x 



Therefore 



h((I + \hf\x'))- 1 + \hf{x'))- 1 ) f(x') 
2h{l- l -h 2 f'{x'f)~ 1 f{x'). 



C Cx , X j X ) C ( X ; X , X ) — C (x , X j X ) C i^X , X , X ) 

— O [X , X , X ' X ) 

= l -VH{x') T {x"-x) 

= 1 -WH(x') T 2h(l~h'f(xry 1 f(x / ) 

= Z -hVH{x') T {Y.^nx')/2) 2n )f{x') 

n=0 

= ^E^')' [{KH"(x'j) 2n K] VH(o 

n=0 

= 0. 



(11) 
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because each matrix in square brackets is antisymmetric. (Each term is an elementary 
Hamiltonian corresponding to a superfluous tall tree.) The expression is a rational function 
of x and h so if its Taylor series in h is zero, the function is zero. Therefore the Kahan method 
has a first integral C(x, x, x'). This can also be written in the symmetric form 

C(x, x, x') = (C(x, x, x') + C(x, x', x'))/2 = C(x, (x + x')/2, x') 

or explicitly as a function of x as 

C(x, x, x') = C(x, x, x + h(l- \hf{x)) ^f(x)) 

= H(x) + l -hVH(x) T (l - hf'(x)y 1 f(x) 
= H(x). 

As in (11), V H (x) T f (x) 2n f \x) = for all n, so we can also write H(x) in a form manifestly 
even in h, 

H(x) = H(x) + h 2 VH(xf(l - h 2 f(x) 2 y 1 f(x)f(x). 

Now consider the case that H{x) is cubic but not homogeneous. We extend it to a 
homogeneous function H(xq, x\, . . . , x n ) so that H(l,x%, . . . ,x n ) = H (x\, . . . , x n ), and extend 
K to K by adding a zero initial row and column, so that i = 0. The linear integral xq is 
conserved by Kahan's method, and KVH(l,x) = KWH(x), so the modified Hamiltonian of 
Kahan's method for KVH(x) reduces to a modified Hamilton of Kahan's method for KWH(x) 



given by the same formula ( 10 1 as in the homogeneous case. □ 



The Kahan map and its conserved quantity H are rational functions of x whose degrees are 
described in the following proposition. When K has full rank, for n — 2 (resp. 3,4), H is degree 
3 over degree 2 (resp. degree 5/2, degree 5/4) and the Kahan map is degree 2/2 (resp. degree 
3/2, degree 4/4). In the planar case, Kahan's method gives a rational map with cubic invariant 
curves. We conjecture that the dynamics of the Kahan map in two dimensions is related to the 
abelian group structure of elliptic curves (and in higher dimensions, to that of abelian varieties) 
as is the case for planar QRT maps [3]. 

Proposition 4. Let H be a cubic in K™ and let K be a rank k antisymmetric n x n matrix. 

(i) The degree of the denominator of H is at most k and the degree of the numerator of H is 
at most k + 3. When k = n the degree of the numerator of H is at most k + 1. 

(ii) The degree of the denominator of the Kahan map is at most k and the degree of the numerator 
is at most k + 1. When k — n the degree of the numerator is at most k. 

Proof, (i) Because the method is linearly covariant we can assume without loss of generality 
that K is in its Darboux normal form 
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where L is k x k. Numbering the blocks of K 1 and 2, the denominator of H is equal to 
det(J - ^/') = det(J - )^hKH"{x)) 



2 v° °/ v^M*) h 22 {x) 

= dei ^I-\hLH ll (x) -\hLH l2 (x) 

= det(7- i/i£/fn(V)). 

Each entry of the matrix LHn(x) is linear in x so the final determinant has degree at most 
k. Next, we write 

~ Hjx) detij - \hf{x)) + \h(VH{x)) T adj(/ - \hf'{x))f{x) 

(X) det(J- \hf'{x)) 

where adj(A) = A^&etA is the adjoint of A. The first term in the numerator has 
degree at most k + 3, VH and / have degree at most 2, and adj(7 — ^hf'(x))f(x) = 
adj(7 — ^hLH\i(x))LH\{x) where each entry in adj(I — \hj'(x)) is the determinant of a 
(k — 1) x (k—1) matrix whose entries are linear in x. Hence the degree of the numerator of 
H(x) is at most k + 3. 

Finally we consider the case k — n. Since k is even, n must be even. First consider the case 
that H is a homogeneous cubic. Then VH(x) — \H"{x)x and H{x) = ^x T H" (x)x. Thus 
in this case we have 

6H(x) = 6H(x) + 2h\7H{x) T (l - ^hf'(xfj f(x) 

= x T (H"{x) + H"{x)(2h- 1 I - KH"(x))- 1 KH"{x)) x 
= x t (h"{x){I- ^hKH"(x))- i y. 

Expanding the matrix inverse using Cramer's rule now shows that the degree of the 
numerator is at most k + 2. The terms of degree k + 2 come from terms in the minors 
of I — ^hKH"(x) of degree k — 1 in x. Every Xi in the matrix I — ^hKH"(x) is multiplied 

by h, thus these terms also have degree k — 1 in h. However, H(x) is an even function of h 
and so these terms must sum to zero. Thus the degree of the numerator is at most k + 1. 
When H (x) is a nonhomogeneous cubic, the terms of degree k + 3 in the numerator of H{x) 
come from the cubic terms in H(x) only, and thus vanish as in the homogeneous case. The 
terms of degree k + 2 are odd in h and hence vanish as before. 

(ii) The proof for the general case follows as above. For the case k — n, we first consider the 
case that H is a homogeneous cubic. Then f(x) — ■^f'(x)x and the Kahan map can be 
written 

x' = x+(l- ^f'(x)y 1 f(x) = (j - \f{x))~\. 

Expanding the matrix inverse using Cramer's rule now shows that the degree of the 
numerator is at most k + 1. In the nonhomogeneous case, the terms of degree k + 2 in 
the numerator come from the cubic terms in H only, and thus vanish as in the homogeneous 
case. 

□ 



Geometric properties of Kahan's method 7 

Examples suggest that there are no other values of n or k other than n = k which lead to a 
reduction in degree. 

Proposition 5. Kahan's method preserves the measure 

dxi A ... A dx n 
dct(7- \hf'{x)) 

for all cubic Hamiltonians on symplectic vector spaces and on Poisson vector spaces with constant 
Poisson structure. 

Proof. Let A = be the Jacobian of the Kahan mapping. Differentiating the mapping |9| 
gives 

^ = af'(x) + 1(1 - 2a)f((x + x')/2)(I + A) + af(x')A 

= af(x) + i(l - 2a)(f(x) + f(x'))(I + A) + af'(x')A 
because f'(x) is linear in x. Solving for A, 

A = (7 - (i + \a)hf\x') -{\- ^hfix))- 1 (/+(! + \a)hf'{x) + (\- \a)hf(x')) . 

Now f'(x) = KH"(x) =: KS where S is symmetric. From Sylvester's determinant theorem, 
det (J + KS) = det (J + SK) = det ((J + SK) T ) = dct(7 - KS). The sum of such matrices has 
the same property, so 

det A = det(/ - (I + \a)hf'{x) -{\- \a)hf'{x')) 
det(J - (i + ia)^/'^') - (| - |o)/»/'(a:)) ' 

This yields invariant measures m(x)dx\ A ... A dx n in 3 cases: 

(i) when a = —1/2 (Kahan's method), m(x) = 1/ det (I — ^hf'(x)); 

(ii) when a = (midpoint rule), m(x) = 1; 

(iii) when a = 1/2 (trapezoidal rule), m(x) = det(I — \hf'{x)). 

□ 

By integrable symplectic map we adopt the definition of Bruschi et al. pQ: a symplectic 
map on a 2n-dimensional symplectic manifold is integrable if it has n functionally independent 
integrals in involution. We will say that leaf-preserving Poisson maps are integrable if the the 
map is integrable on each leaf. 

Corollary 6. Kahan's method yields an integrable mapping of the plane when applied to any 
canonical Hamiltonian system in the plane with cubic Hamiltonian. Kahan's method yields an 
integrable mapping of M 3 when applied to any Poisson system on R 3 with constant Poisson 
structure and any cubic Hamiltonian. 

Proof. A measure and a first integral are sufficient for integrability in the plane. The odd- 
dimensional case with constant K has a linear Casimir which is conserved by the method, 
reducing the situation in this case to two dimensions on each level set of the Casimir. □ 

Corollary 7. When n — 2 and H is a homogeneous cubic, H{x)=H{x)/det{I-\hf(x)) and 
Kahan's method preserves the h-independent measure {dx\ A dx2)/H(x). 
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Proof. We have 

6det(J- ~hf'(x))H(x) = det(J - ~hf (x))x T (h" (x)(I - ^hKH"[x))- x ^x 
= x T (^H" {x) &d]{I - -hKH"(x)fjx 
= x t (h"(x)(I- adj(~ hKH"(x)))x (because n = 2) 
= x T (H"(x) - ^hH"{x){KH"(x))~ 1 det(KH"(x))^x 

= x t (h"(x) - -hKdet(KH"(x)f)x 

= x T H"{x)x (because K T = -K) 
= 6H(x). 

Any map that preserves a measure n(x) and an integral I(x) also preserves the measure I(x)n(x). 
Taking fi = dx\ A dx^j det(7 — \hf{xj) and I(x) — 1/H(x) gives the /i-independent measure 
(dxi Adx 2 )/H(x). " □ 



4. Discussion 

Level sets of H are shown in Figure 1 for H — \{q 2 +p 2 ) + q 2 p— \p 3 - Notice that the separatrices 
persist (and are unchanged) for alln, but that the singular set q 2 + p 2 = | + moves in from 
infinity as h increases and alters the topology of the level sets. For h < y4/3 the topology of 
the bounded orbits is unaltered. 

The bounded orbits of Figure [T] are symmetric, so for the following numerical experiments 
we used H = p — p 3 + q 2 — q 3 , which has bounded, nonsymmetric orbits, an elliptic fixed point 
at (q,p) — (2/3, and a separatrix meeting (q,p) — (0, l/y/3). Level sets of H for this case 

are shown in Figure [2] 

Numerical experiments strongly indicate that the following observations hold. 

(i) No other method of the family fcty has a modified Hamiltonian when H is cubic, apart from 
the known cases a = 0, ±|, and~| (see Figure [3]). 

(ii) The midpoint and trapezoidal rules do not have a first integral for all cubic H (even though 
they do have a formal invariant close to H ) (see Figure [4]) . 

(iii) Simpson's method is not measure- preserving for all cubic H. (For H = p — p 3 + q 2 — q 3 , 
a numerical calculation finds eigenvalues 1, A of periodic orbits, with A =/= 1, contradicting 
measure preservation.) 

(iv) Kahan's method does not preserve any symplectic form in dimension > 4 for all cubic H. 
(A numerical calculation of periodic points finds eigenvalues that do not occur in A, 1/A 
pairs. Proposition [2] establishes this for a limited class of symplectic forms.) 

(v) Compositions of Kahan's method with different step sizes do not have a modified 
Hamiltonian when H is cubic (see Figure [5]) . 

Our results are significant and novel for the study of both the integrability of the mappings 
produced by Kahan's method and for the study of the geometric properties of Runge-Kutta 
methods: 
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Figure 1. Top left: Level sets of H = ~ (<j 2 + p 2 ) + q 2 p — ^p 3 (the so-called Henon— Heiles 

potential). Same level sets of the conserved quantity H of Kahan's method for h = 1/3 (top 
right); h = 2/3 (bottom left) (the jagged circle q 2 + p 2 = j + A = 1.58 indicates H = oo, 
on which initial conditions are mapped to infinity — for h = 1/3 the circle has radius 3.04 and 
is out of view); and h — > oo (bottom right). Note that Kahan's method preserves the 3- fold 
discrete symmetry of H, because as a Runge-Kutta method it preserves all affine symmetries. 





Figure 3. Measured rate of energy drift for H = p — p 3 + q 2 — g 3 for Runge— Kutta methods 
x' = x + h(af(x) + (1 — 2a)f((x + x')/2) + a/(x')), varying the parameter a. The step size is 
h = 0.3 and the initial condition is q = 0.323, p = l/\/3. All methods have an approximate 
modified energy up to h 4 . The energy drift is measured by fitting a straight line to this 
modified energy over 2 X 10 6 time steps. Only the four methods identified by the analysis 
(a = —1/2, Kahan; a = 0, midpoint; a = 1/6, Simpson; and a = 1/2, trapezoidal) show no 
energy drift by this measure. 




Figure 4. Portion of the phase portrait of the midpoint rule with step size h = 0.3 applied 
to H = p — p 3 + q 2 — q 3 . The observed chaotic bands and island chains indicate that it does 
not have a conserved quantity. 
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Figure 5. Suzuki's 3-stage, 4th-order composition applied to Kahan's method shows a 
comparatively rapid energy drift, indicating that there is no conserved quantity. Here 



H - 



-p-p" 



h = 0.2, and (<?o,Po) = (0.323, 1/V3). 



• First, our results explain the intcgrability of the map obtained when Kahan's method is 
applied to some of the examples of [ID] : their Eq. (4.2) (H = y 2 /2 — 2x 3 + ax); Eq. (5.4) 
(H = y(3x 2 — y)); Eq. (8.1) (Volterra chain in M 3 , H = xix 2 x 3 , constant K, integral H 2 
in Eq. (8.6) is a function of our H and the Casimir xi + x 2 + x$); Eq. (9.1) (Dressing 
chain in K 3 , H = (x\ + x 2 ){x2 + £3X^3 + £1) — J2i a i x i, constant K). Our results explain 
the invariant measure and cubic integral of their Eq. (11-1) (three wave system in C 3 , 
H = Z1Z2Z3 + Z1Z2Z3), the invariant measure for the family of systems in their Prop. 1, and 
the linear integrals throughout [TO] . 

• Second, our results (e.g. Corollary[6]) systematically produce new integrable cases of Kahan's 
method. 

• Third, we have shown that Kahan's method in dimension 4 and greater provides examples of 
maps with nonlinear integrals and conserved measures unrelated (in general) to intcgrability 
or obvious symmetries, again a novel feature. For example, our results imply that Kahan's 
application of the method to the Korteweg-de Vries equation in [5] preserves a measure and 
a modified energy (but the higher order compositions of the method in [5] probably do not). 

• Fourth, we have shown that Kahan's method has novel properties previously unknown 
amongst Runge-Kutta methods, indeed amongst all B-series. It is known that B-series 
methods cannot conserve the measure dx\A. . .Adx n even for linear vector fields 0; Kahan's 
method circumvents this by conserving a modified measure. It is a novel conjugate-to-energy 
preserving method for cubic H . In the plane, it is also conjugate to symplectic. Thus, while 
no conjugate to symplectic methods are known that are also energy preserving in general, 
here we have one that preserves at least a modified energy, and preserves it exactly (not 
merely as a formal invariant). 

On the other hand, there are open questions in all of these areas. While it was already 
suggested in |10j that there could be an underlying 'integr ability mechanism' unifying the 
integrable cases, here we have unified only some of these. In addition the hoped-for unification 
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should now be extended to include non-integrable cases preserving a measure and/or some 
integrals as well. On the numerical side, it is not known precisely which Runge-Kutta or B- 
series methods share the properties of Kahan's method, if any are higher order integrators, or if 
any are conservative for nonquadratic (e.g. other polynomial) vector fields. 
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